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ABSTRACT 

Potential magnetic field solutions can be obtained based on the synoptic magnetograms of 
the Sun. Traditionally, a spherical harmonics decomposition of the magnetogram is used to 
construct the current and divergence free magnetic field solution. This method works reasonably 
well when the order of spherical harmonics is limited to be small relative to the resolution of the 
magnetogram, although some artifacts, such as ringing, can arise around sharp features. When 
the number of spherical harmonics is increased, however, using the raw magnetogram data given 
on a grid that is uniform in the sine of the latitude coordinate can result in inaccurate and 
unreliable results, especially in the polar regions close to the Sun. 

We discuss here two approaches that can mitigate or completely avoid these problems: i) 
Remeshing the magnetogram onto a grid with uniform resolution in latitude, and limiting the 
highest order of the spherical harmonics to the anti- alias limit; ii) Using an iterative finite differ- 
ence algorithm to solve for the potential field. The naive and the improved numerical solutions 
are compared for actual magnetograms, and the differences are found to be rather dramatic. 

We made our new Finite Difference Iterative Potential-field Solver (FDIPS) a publically avail- 
able code, so that other researchers can also use it as an alternative to the spherical harmonics 
approach. 

Subject headings: Sun, magnetogram, potential field, spherical harmonics, finite-difference method 



1. Introduction 

Magnetograms provide the radial magnetic field on the visible surface of the Sun. The actual mea- 
surement is for the line-of-sight component of the magnetic field, which is then transformed into the radial 
component assuming an (approximately) radial field near the solar surface. As the Sun rotates, the individ- 
ual magnetograms can be combined into a synoptic magnetogram that covers the whole spherical surface. 
Synoptic magnetograms are provided by many observatories, including Wilcox Solar Observatory (WSO), 
the Michelson Doppler Imager (MDI) instrument on the Solar and Heliospheric Observatory (SOHO), the 
Global Oscillation Network Group (GONG), Solar Dynamic Observatory (SDO) and the Synoptic Opti- 
cal Long-term Investigations of the Sun (SOLIS) observatory. Today's magnetograms contain hundreds to 
thousands of pixels along each coordinate direction. These magnetograms can be used to extrapolate the 
magnetic field into the solar corona. 

The simplest model ( Schatten. Wilcox, fc Ness 19691) assumes a current- free, in other words potential, 



magnetic field that matches the radial field of the magnetogram on the surface, while it satisfies a simple 
boundary condition at the outer boundary at some radial distance R. The outer boundary condition is 
usually taken at R = 2.5 (solar radii), and a purely radial field is assumed at this "source surface". 
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Mathematically the problem is the following: given the magnetogram data that defines the radial component 
of the magnetic field as M{6, (f>) at r = 1 Rs, find the scalar potential $ so that 

V-(V$) = (1) 



dr 



= M{0,(b) (2) 

K^R = (3) 



Here 9 e [0, tt] and <j) G [0, 27r] are the co-latitude and longitude coordinates, repectively. Once the solution 
is found, the potential field solution is obtained as 

B = V$ (4) 

and it will trivially satisfy both the divergence-free and the current-free properties 

V • B = V • (V$) = (5) 
V X B = V X (V$) = (6) 

We note that the current is only zero inside the domain. If the solution is continued out to r > R with a 
purely radial magnetic field Br{r > R) = {R/r)^Br{R), there will be a finite current at r > i?, on the other 
hand, the divergence will be zero for all r > 1. 



The potential field solution is often obtained with a spherical harmonics expansion (|Altschuler et al 



19771 ). Here we briefly summarize the procedure in its simplest possible form. The base functions ipnm are 



the spherical harmonic functions Ynm multiplied with an appropriate linear combination of the corresponding 
radial functions r" and r~"~^ so that the boundary condition (pnm{R, = is satisfied: 

^n^{r,e,(b) - (r" - i?2"+ir-"-i)y„,„(0,0) (7) 

The indexes n > and m (|m| < n) are the integer degree and order of the spherical harmonic function, 
respectively. The (pnm functions are solutions of the Laplace equation ([Ij , satisfy the boundary condition at 
r — R, and they form an orthogonal base in the 9, (p coordinates. The magnetic potential solution can be 
approximated as a linear combination of the base functions 

N n 

$(r-,61,0) = ^ fnm^nmir,9,(j)) (8) 

n—1 m—~n 

where N is the highest degree considered in the expansion and the n — m — harmonics is not included, as 
it corresponds to the monopole term. The coefficients fnm can be determined by taking the radial derivative 
of equation ([8]) and equating it with the magnetogram radial field at r = 1 : 



N n 



N 



n—1 m— — n r— 1 n—1 m— — n 

Exploiting the orthogonality of the base functions, we can take the inner product with Ynm to determine 

= A I ^( \uwn+u r sm9 d<j,M{9, (f>)Y„„,{9, 0) (10) 

where the l/(47r) coefficient results from the normalization of the spherical harmonics. An alternative 
approach of obtaining the harmonic coefficients /„,„ is to employ a (least-squares) fitting procedure in 
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equation ([9]). This is much more expensive than evaluating the integral in (jlOp . but it can be more robust 
if the magnetogram does not cover (well) the whole surface of the Sun. 

Using the spherical harmonic coefficients the potential can be determined on an arbitrary grid using 
and the magnetic field can be obtained with finite differences. Alternatively, one can calculate the gradient 
of the base functions analytically and obtain the magnetic field as 

N n 

B(r,0,0) = ^ /nmV(^„™(r,0,0) (11) 

n— 1 ni—~n 

for 1 < r < R. Spherical harmonics provide a computationally efficient and very elegant way of solving the 
Laplace equation on a spherical shell. However, one needs to be cautious of how the integral in equation 
(|10P is evaluated, especially when a large number of harmonics are used in the series expansion. 

We will use the GONG synoptic magnetogram for Carrington Rotation 2077 (CR2077, from November 
20 to December 17, 2008) as an example to demonstrate the problem. The magnetogram contains the radial 
field on a 180 x 360 latitude-longitude grid on the solar surface. The grid spacing is uniform in cos 9 (or sine 
of the latitude) and in longitude (j). Figure [1] shows the radial field. 

Section 2 discusses the naive and more sophisticated ways of obtaining the potential field solution with 
spherical harmonics. Section 3 describes an alternative approach using an iterative finite difference. The 
various methods are compared in the final Section 4, where we also demonstrate the ringing effect that can 
arise in the spherical harmonics solution, and we draw our conclusions. 



2. Potential Field Solution Based on Spherical Harmonics 

To turn the analytic prescription given in the introduction into a scheme that works with real magne- 
tograms, one has to pick the maximum degree N, and evaluate the integrals in equation (jlOp for each pair 
of n and m up to the highest order. The resulting /„„ coefficients can be used to construct the 3D potential 
magnetic field solution at any given point using equation (jlip . 

2.1. Naive Spherical Harmonics Approach 

The simplest approximation to equation (jlOp is a discrete integral using the original magnetogram data: 

^ Ng W« 

= 4^[n+ (n + l)fl^"+i1 ^^^^'"'^^'^^'^^^"^"^•^"'"^^"^^•^ ^^^^ 

where the Mij is the radial field in a pixel of the Ng by iV^ sized magnetogram. The pixel is centered at 
the {9i,(f)j) coordinates, and the area of the pixel is given by (A cos 6')i(A(?!))j . 

Unfortunately, the uniform cos 9 mesh used by most of the magnetograms is not at all optimal to 
evaluate the integral in equation (jlOp . In fact this procedure will only work with maximum order N that is 
much less than Nq. Figure [2] shows the Pgo.o associated Legendre polynomial discretized in different ways. 
The red curve shows the discretization on 180 grid points that are uniform in cos 9. Clearly the red curve 
is a very poor representation near the poles, where cos9 = ±1. This is important, because the amplitude 
of the Legendre polynomial is actually largest near the poles. This means that the orthonormal property 
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is not satisfied in the discrete sense, and the coefRcients obtained with equation (|12p are very inaccurate. 
The Legendre polynomial can be represented much better on a uniform grid (shown by the green curve in 
Figure [5]), as we will discuss below. 

A clear signal of this problem is that the amplitudes of the higher order spherical harmonics are not 
getting smaller with increasing indexes n and to, i.e., the harmonic expansion is not converging. The black 
line in Figure [3] shows the amplitudes Sn = X)™ fnm which is oscillating wildly for n > 60 for this 360 x 180 
magnetogram. The oscillations are almost exclusively due to the f„n coefficients, the m coefficients are 



well behaved. This plot can be directly compared with Figure 15 in lAltschuler et al.l ([1977|), where the power 
spectrum is more-or-less exponentially decaying. We believe that the reason is that these authors used a 
least-square fitting to the line-of-sight (LOS) magnetic field instead of calculating the spherical harmonics 
from the radial field as shown above. While the two methods are identical analytically (assuming that the 
LOS and radial fields correspond to the same solution), the use of least-square fitting mitigates the lack of 
orthogonality among the discretized Legendre funtions, while the naive approach described above heavily 
relies on the orthogonality property. 

Given the non-converging series expansion, the resulting potential field will be very inaccurate in the 
polar regions and will have essentially random values depending on the number of spherical harmonics used. 
This is demonstrated by Figure |3] that shows the radial magnetic field reconstructed with various number 
of harmonics using the original magnetogram grid. One would expect the radial component of the potential 
field to reproduce the magnetogram shown in Figure [TJ Instead, we find that the solution deviates strongly 
in the polar region if the harmonics expansion is continued above N = 60. For TV = 60 (or lower) the solution 
looks reasonable, but strongly smoothed due to the insufficient number of harmonics. This is most obvious 
around the active regions in the top panel of Figure SI 

We note that these numerical errors are not related or comparable to the observational uncertainties of 
the magnetograms, which are usually also quite large in the polar regions. The observational uncertainties 
are essentially unavoidable but are within some well-understood range. On the other hand, these numerical 
artifacts are definitely avoidable while the errors are essentialy unbounded if one uses too many harmonics. 



2.2. Spherical Harmonics with Remeshed Magnetogram 

One can get much more accurate results if the magnetogram is remeshed to a grid that is uniform 
in the co-latitude 9, has an odd number of nodes, and contains th e two poles ^ = and 9 = tt. In fact 



this is the standard grid used in spherical harmonics transforms (e.g. ISuda fc Takamil (|2001[ )) and it is often 
referred to as using the Chebyshev nodes, since the uniform 9 grid points correspond to the Chebyshev nodes 
in the original cos 9 coordinate which is the argument of the Legendre polynomials. Figure [2] shows that 
the Legendre polynomial is much better represented on the uniform 9 grid than on the uniform cos 9 grid. 
Remeshing the magnetogram introduces some new adjustable parameters into the procedure: the number of 
grid cells Ng on the new mesh, and the interpolation procedure. 

If the remeshing is done with the same number of grid points as is in the original magnetogram grid, 
the latitudinal cell size at the equator will be a factor of tt/2 larger than in the uniform-cos grid. On 
the other hand, the uniform-^ grid will contain many more points than the original in the polar regions, 
so the interpolation procedure may create some unwanted artifacts. To maintain the resolution of the 
original data around the equator, we set Ng to {'k/2)Ns rounded to an odd integer. For the remeshing 
we chose a simple linear interpolation procedure, and it works satisfactorily, but one could certainly use 
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higher order interpolation procedures, such as splines. Before doing the interpolation, we add extra grid cells 
corresponding to the north and south poles of the magnetogram grid, and the values at these two extra cells 
are set as the average of the pixels around the poles: 



Mn = 



^Em.. (13) 

1 = (14) 



The co-latitude coordinates of the uniform-0 mesh arc 



« - 1 

e' = TT— (15) 

for i' — 1 . . . Ng. We use simple linear interpolation from the extended magnetogram mesh to the uniform 
mesh: 

M;, = aM,,j- + (1 - a)M,-+ij (16) 
where the i index is determined so that 9i < 6[, < Oi^i and 

Q' 



(17) 



Oi+l — Oi 

Finally the spherical harmonics coefhcients are determined with the integral approximated as 

" 47r[n+(n+l)j;2n+i1 J2 J2 e.u;^(A0)jM^y„™(0^, 0^) (18) 

where ei = e^r' = 1/2 and = 1 for all other ind exes. The Wi coefficients are the Clenshaw-Curtis weights 
Icienshaw fc Curtij|l96ci[ IPotts. Steidl. fc TaschaJQQa ) defined as 

"'^-;^E^'^ifc^'=°^(2^^^') (19) 

where H — {Ng — l)/2 and eg = = 1/2 and ej. = 1 for all other indexes. We note that for Ng order of 10 
or more 

|cos6'^^ -cos6'^ 

« J ^ (20) 

where 1+ — min(i + 1, Ng) and i- — max(i — 1, 1) are the indexes of the neighboring cells, or the cell itself 
at the poles. 

Using the proper grid allows us to use larger nu mber of spherical har monics. N is limited only by 
the TV < 2iV^/3 and TV < N^/S ahas-free conditions |Suda fc TakamilboOll ). For our example 180 x 360 
magnetogram, we remesh it a to 283 x 360 uniform-^ grid, and we can obtain accurate solutions up to 
N ~ 120 degree spherical harmonics. Figure [5] shows the potential field solution obtained with the remeshed 
magnetogram grid with N = 120. Compared with the naive method, the solution is much more reasonable 
in the polar regions. There is some smoothing when compared to the magnetogram shown in Figure [U most 
obvious near the active regions. 

While the remeshing is definitely a big improvement over using the original magnetogram, it would be 
nice to be able to use the original magnetogram data without remeshing and interpolation. The next section 
shows that this can be easily achieved with a finite difference solver. 
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3. Finite DifTerence Iterative Potential Field Solver (FDIPS) 

The Laplace equation ([Ij with the boundary conditions ^ and ([3]) can be solved quite easily with an 
iterative finite difference method. The advantage of finite differences compared with spherical harmonics 
is that the boundary data given by the magnetogram directly effects the solution only locally, while the 
spherical harmonics are global functions, and their amplitudes depend on all of the magnetogram data. If 
the magnetogram contains large discontinuities, we expect the finite difference scheme to be better behaved. 

The finite difference method has advantages if the solution is to be used in a finite difference code on the 
same grid, because one can guarantee zero divergence and curl for the magnetic field in the finite difference 
sense. The solution obtained with the spherical harmonics has zero divergence and curl analytically, but not 
on the finite difference grid, which may severly underresolve the high order spherical harmonic functions in 
some regions (see Figure [2]). 

The finite difference m ethod was applied to the solar potential magnetic field problem as early as 1976 



([Adams fc PneumanI (|1976I )). but the method was limited by the computational resources available at the 
time. Solving a 3D Laplace equation on today's computers is an almost trivial problem. We implemented 
the new Finite Difference Iterative Potential-field Solver (FDIPS) code in Fortran 90. The serial version 
does not require any external libraries, while the parallel version uses the Message Passing Interface (MPI) 
library for communication. FDIPS can solve the Laplace equation on a 150 x 180 x 360 spherical grid to 
high accuracy on a single processor in less than an hour. The parallel code can solve the same problem in 
less than 5 minutes on 16 processors. 

We briefiy describe the algorithm in FDIPS. We use a staggered spherical grid: the magnetic field is 
discretized on cell faces while the potential is discretized at the cell centers. We use one layer of ghost cells to 
apply the boundary conditions so the cell centers are located at , 6j , 0^ with i = O...Nr + l,j = O...Ne + l 
and k = . . . + 1. The 9j and 0^ coordinates of the real cells are given by the magnetogram, while the 
ghost cell coordinates are given by 6*0 = —6*1, ^^jVe+i = 27r — ^at^, (/)o — 4>n^, and 4>n^+i = ipi- We allow for a 
non-uniform radial grid extending from r = 1 to i?, but for sake of simplicity in this paper a uniform radial 
grid is used with = I + {i - 1/2) Ar with Ar = {R - r)/Nr. 

The radial magnetic field components are located at the radial cell interfaces at (?'i+i/2, ^j, '/'fc) where 
'''i+i/2 = (^i + ^i+i)/2 for i = . . . Nr, j = 1 . . . Ng and k ^ 1 . . . N^. Similarly the latitudinal components 
are at r^, 9j^i/2, 4>k with cos6'j_|_i/2 (cos 6j + cos6'j+i)/2, and j — . . . Ng. Note that the interface is taken 
half-way in the cos^ coordinate and not in 9, because this makes the cells equal area when the magnetogram 
is given on uniform cos 9 grid. Finally the longitudinal field components are located at (r^, 9j, 4'k+i/2) where 
'/>fc+i/2 = (0fe + 0fe+i)/2 for fc = . . . A^0. 

The staggered discretization keeps the stencil of the Laplace operator compact and it makes the boundary 
condtions relatively simple. The magnetic field is obtained as a discrete gradient of $: 



Br,i+l/2,j,k 



Be,i.j + l/2,k 



Ar 

sin6'j+i/2($i,j+i,fe - $i,j,fe) 



r^Acos^ 
r,; sin 9jA(f) 



„ _ ^i,j.k+i - ^i.j,k 



Note the sin 0/ A cos 6* factor in the 9 derivative for the uniform cos 9 grid. For uniform 9 grid this is replaced 
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with 1/A6'. The divergence of the magnetic field, i.e. the Laplace of <i> is obtained as 

„ /VT2^^ 7'?+l/2-Sr,i+l/2J,fe - '^i-l/2,j.k^r,t-l/2,3,k 
= ( V ^)^.J,k = ^3^^^ 

^ sin6'j+i/2i?6i,i,j+i/2,fe - sin6'j_i/2-B6),ij_i/2,fc 
A cos 6* 

-B0,i,j,fc+i/2 - -S0,ij,fc-i/2 .22^ 
r,;sin6ljA0 ^ ^ 

Again 1 /A cos 9 is used for the uniform cos 9 grid, while on the uniform 9 grid this is replaced with l/(sin 9A9). 
The magnetogram boundary condition is applied by setting the inner ghost cell as 

$Oj,fc = *ij,fc-ArMjfe (23) 

where Mj = Mj,k — M is the magnetogram with the average field (i.e. the monopole due to observation 
errors) 



p^(Acos0)j(A0)feM,-fc (24) 



An 
j.k 

removed. The zero potential at rN^+1/2 = R is enforced by setting the ghost cell as 

^N^ + l,j,k = -^N,.J,k (25) 

The boundary conditions at the poles are a bit tricky. Cells (i,l,fc) and {i,l,k') are on opposite sides of 
the north pole if k' = mod(A: — 1 + N^/2,Nij,) + 1. Therefore the ghost cells in the 9 direction are set as 
'^i.o,k = $i,i.fc' and ^i.Ne+i,k = ^i,Ng.k' ■ We note here that iV^ is assumed to be an even number. The 
periodic boundaries in the cf) direction arc simple: ^i,j.o — ^ij.N^ and ^ij^N^^i — 



We need to find ^ij^k that satisfies the discrete Laplace equation p2p with the boundary conditions 
applied via the ghost cells. The initial guess is $ = which provides a non-zero residual because of the 
inhomogeneous boundary condition at the inner boundary applied by equation (|23p. We use this residual 
with a negative sign as the righ-hand-side of the Poisson equation (V^<I>)ij^fe — Rij.k, and use *&oj,/c = ^i.j,k 
as the new homoge neous inner b oundary condition instead of (|23p . We use the Krylov-type iterative method 
BiCGSTAB ( van der Vorstlll992h to find the solution. The linear system is preconditioned with an Incomplete 



Lower-Upper decomposition (ILU) preconditioner to speed up the convergence. We use ILU(O) with no fill- 
in compared to the original matrix structure, so the preconditioner is a diagonal matrix, but its elements 
depend on all elements of the original matrix. We have implemented a serial as well as a parallel version of 
the algorithm. In the parallel version the preconditioner is applied separately in each subdomain. FDIPS 
finds an accurate (down to 10^^° relative error) solution on a ISO'^ grid in less than 1000 iterations. Even 
running serially, this takes less than an hour on today's computers. 

Once the solution is found in terms of the discrete potential ^i,j,k, we apply the original boundary 
conditions including ((23|) and calculate the magnetic field with equation (|2T|) . The divergence of the magnetic 
field will be zero to the accuracy of the Poisson solver. The curl of the magnetic field will be zero in a finite 
difference sense simply because it is constructed as the discrete gradient of the potential. The boundary 
condition at the inner boundary is also satisfied exactly: -Br,i=i/2j".fe = -^j k- Averaging rB^ and rBg to the 
— R position at the outer boundary also gives exactly zero tangential fields due to the equations ()25|) 
and ([2T|l . Depending on the application, we may interpolate the potential magnetic field onto a collocated 
grid, or use it on the original staggered grid. 
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Figure [6] shows the solution of the magnetic field obtained with the finite difference solver FDIPS on 
a 150 X 360 x 180 grid. Since we use the same uniform-cos 6 grid as the magnetogram, the obtained radial 
field is identical with the magnetogram at r = 1. The tangential components agree well with the remeshed 
spherical harmonics solution shown in Figure [SJ It took 1166 iterations to get a solution with a relative 
accuracy of 10^^". The run time was almost exactly one hour on a 2.66 GHz Intel CPU. 



4. Discussion and Conclusions 

We have discussed various ways to obtain the potential field solution based on solar magnetograms. 
While spherical harmonics provide an efficient and elegant method, there are some subtle restrictions that 
require attention. If one wants to use many spherical harmonics (the same order as the number of magne- 
togram pixels in the colatitude direction), the magnetogram data on the Ng x iV^ grid has to be remeshed 
onto a uniform-^ grid with Ng x points, Ng must be an odd number, and the new grid must include both 
poles. After the remeshing the maximum degree of harmonics N is only limited by the anti-alias limit to 
min(2A^g/3, N^/3). We used a simple linear interpolation for the remeshing. 

The remeshing can be avoided by the use of a 3D finite difference scheme. One can use the original 
magnetogram grid, and the only freedom is in choosing the radial discretization. The finite difference scheme 
provides a solution that is fully compatible with the boundary conditions, and the solution has zero divergence 
and curl in the finite difference sense. 

Figure [7] compares the solutions obtained with the three methods along the radial direction for a fixed 
latitude 9 = 76.5° and longitude (p = 30°. The spherical harmonics series were truncated at = 120 for both 
the naive and remeshed methods. The naive spherical harmonics algorithm gives incorrect results close to 
the solar surface where the high order harmonics dominate. This is most obvious for the radial component, 
which is given at r = 1 by the magnetogram, and it is exactly reproduced by the finite difference scheme. 
The latitudinal component at r = 1 is also very different from the values given by the remeshed harmonics 
and the finite differences. The latter two methods agree reasonably well with each other. For radial distances 
above r = 1.05 all three methods agree quite well. 

So far we restricted our example to a GONG magnetogram taken at the solar minimum. If one uses an 
MDI magnetogram during solar maximum, the largest magnetic fields are much stronger (order of 1000 G) 
and the resolution of the magnetogram is much finer (order of 1000 pixels). The finer magnetogram resolution 
allows going to larger number of harmonics, even when using the original magnetogram grid (naive approach) . 
But the strong and sharp gradients in the magnetogram will bring out another problem with the spherical 
harmonics approach, the ringing effect. The ringing is due to the so-called Gibbs phenomenon: the step- 
function like magnetogram data results in high amplitude high o rder h armonics in Fourier space. The ringing 
effect and other artifacts are discussed in great detail by Tran ( 20091 ). 



Figure [S] demonstrates this effect on the 3600 x 1080 resolution MDI magnetogram for Carrington 
Rotation 2029 (from April 21 to May 18, 2005), with the maximum radial field strength around ±3000 G. The 
remeshed harmonics method with A^ = 90 is compared with the finite difference method on a 150 x 360 x 180 
grid (the magnetogram data is coarsened to a 360 x 180 grid). In the spherical harmonics solution the 
ringing is very clearly visible around the active regions, both in the radial and latitudinal components. The 
finite difference scheme, on the other hand, shows no sign of ringing in either components. This is obvious 
for the radial component, which simply coincides with the coarsened magnetogram, but for the latitudinal 
component it is due to the fact that the finite difference solution of the Laplace equation does not suffer from 
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ringing artifacts even for discontinuous boundary data. For the spherical harmonics approach the ringing 
becomes weaker with increased number of harmonics, but it is quite apparent even for N — 180 (not shown). 
The results of the remeshed and naive harmonics methods are essentially the same up to = 180, i.e. the 
ringing is not due to the remeshing of the magnetogram. 

In terms of computational efHciency, a good implementation of the spherical harmonics scheme is much 
faster than the finite difference scheme. In fact, it may be more costly to construct the potential field 
solution on a 3D grid from the spherical harmonics coefficients than obtaining the coefficients themselves. 
Our Fortran 90 code can obtain the spherical coefficients up to = 60, 90 and 120 degrees in 1, 1.8 and 
3.3 seconds, respectively, while the reconstruction of the solution on the 151 x 361 x 180 grid takes 5, 12, 
and 20 minutes, respectively. All timings were done on a single 2.66 GHz Intel CPU. The reconstruction 
cost can be improved by running the code in parallel, and/or truncating the series in parts of the grid where 
the higher order harmonics have a negl igible contribution. We also n ote that going beyond about A^ = 360 
harmonics becomes fairly complicated ( Potts. Steidl. fc Tasche 1998[) . 



The computational cost of the finite difference scheme scales with the number of grid cells and the 
number of iterations. The number of iterations is fairly constant for multigrid type methods, but for the 
Krylov sub-space schemes it grows with the problem size, although slower than linearly. The finite difference 
scheme can be sped up by parallelizing the code, which is fairly straightforward for the Krylov subspace 
schemes. Since we limit the ILU preconditioning to operate independently on the subdomains of each 
processor, the preconditioner becomes less efficient as the number of processors increases, which results in 
an increase in the number of iterations. To minimize this effect, the parallel FDIPS code splits the grid in 
the 9 and (j) directions only, so the subdomains in each processor contain the full radial extent of the grid. 
Our experiments confirmed that using this domain decomposition, the number of iterations indeed does not 
depend much on the number of processors. Our largest test so far involves a 450 x 540 x 1200 grid with 30 
times more cells than the 150 x 180 x 360 grids discussed in most of this paper. For the large problem we 
need about 8,500 iterations to reach the 10"^" relative accuracy, a factor of 9 increase relative to the smaller 
problem. Using 108 CPU-s, the solution is obtained in about 5.3 hours. 

Despite the various limitations, for some applications the spherical harmonics approach may still be 
preferred. For example if the solution is needed to obtain a spherical power spectrum of the solar magnetic 
field. If the solution is to be used in a finite difference code, the finite difference solution is probably 
preferable. We are using the FDIPS code to generate the potential field solution as the initial field for our 
solar corona model ( van der Hoist et al. 2010[) . 



This paper attempts to call the attention of astrophysicists and solar physicists to the limitations and 
potential pitfalls of using the spherical harmonics approach to obtain a potential field solution. The spherical 
harmonics representation of the potential field solutions are available from several synoptic magnetogram 
providers, although the details of the method used to obtain the spherical harmonics is not always clear. A 
spherical harmonics based PFSS package implemented in IDL is available as part of the Solar-Soft library 
(jhttp : //www. Imsal . com/solarsof t). This package uses the magnetogram remeshing technique either onto 
the Chebysliev (uniform-0) or the Legendre collocation points. 

We are not aware of any publically available code that uses finite differences to solve this particular 
problem. To allow other researchers to use and compare the two approaches, we make our finite difference 



code FDIPS publically available at the http : //csem . engin . umich . edu/f dips/| website. 
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Fig. 1. — Synoptic magnetogram for Carrington Rotation 2077 obtained by GONG. The radial magnetic 
field at the photosphere is shown in the range —15 to +15 Gauss to show more detail. The color scale is 
saturated in the active regions where the largest and smallest values are — 45.4G and 27.9 G. 



- 12 - 




-0.3 E , , , , \ , , , , \ , , , , \ , , , , 3 

-1.00 -0.95 -0.90 -0.85 -0.80 

cos Theta 



Fig. 2. — Discrete representations of the Pgco associated Legendre polynomial as a function of cos 6* near 
the "south pole" at cos^^ = — 1. The black curve shows an accurate representation with 1800 grid points 
uniformly distributed in the [—1,1] range. The red curve represents the polynomial on 180 grid points 
uniform in cos0, while the green curve uses 181 grid points that are uniformly distributed in 9 including the 
poles at 9 = and 9 = n. 
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Fig. 3. — Power spectrum Sn — J2m fnm of the spherical harmonics expansion with the original (black line) 
and remeshed (red line) magnetograms. 
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Fig. 4. — Comparison of the radial magnetic field at r = 1 using the spherical harmonics expansion on the 
original magnetogram grid up to = 60 (top), iV = 90 (middle) and N = 120 (bottom) order. These plots 
should be compared with the magnetogram shown in Figure [T] Note that the magnetic field in the polar 
regions is completely wrong for = 90 and TV = 120. The = 60 solution is reasonable, but the fine details 
are smoothed out. 
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Fig. 5. — Magnetic field solution at r = 1 using the remeshed spherical harmonics up to TV = 120 degree. 
The radial compnents (top panel) well reproduces the magnetogram shown in Figure [TJ although some of 
the details are slightly smoothed out. 
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Fig. 6. — Magnetic field solution at r = 1 using the finite difference code FDIPS on a 150 x 360 x 180 grid. 
The radial component (top panel) agrees exactly with the magnetogram shown in Figure [1] except for the 
removal of the average field (the monopole). The tangential components (middle and bottom panels) agree 
well with the remeshed spherical harmonics solution shown in Figure [5j 
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Fig. 7. — The radial and latitudinal components of the magnetic field along the radial coordinate at fixed 
76.5° latitude and 30° longitude. The solutions are obtained with the naive (blue line) and remeshed (red 
line) harmonics approaches with N ~ 120, as well as with finite differences (black line). 
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Fig. 8. — Radial and latitudinal components of the potential field solution at r = 1 for the MDI magnetogram 
for CR2029 obtained with remeshed spherical harmonics algorithm with maximum order TV = 90 (left) and 
with FDIPS using a 150 x 360 x 180 grid (right). The color scale is saturated in the active regions. 



